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Abstract 

We consider a reaction-diffusion system including discontinuous hysteretic relay 
operators in reaction terms. This system is motivated by an epigenetic model that 
describes the evolution of a population of organisms which can switch their phenotype 
in response to changes of the state of the environment. The model exhibits formation 
of patterns in the space of distributions of the phenotypes over the range of admissible 
switching strategies. We propose asymptotic formulas for the pattern and the process 
of its formation. 
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1 Introduction 

In this paper, we consider a reaction-diffusion equation that includes discontinuous two- 
state two-threshold relay operators. This equation was proposed as a model for dynamics 
of a population of bacteria that can switch their phenotype in response to changes of the 
environment H2G2]. The model exhibits formation of patterns in the space of distributions of 
the two phenotypes over an admissible range of the switching thresholds. The convergence 
of solutions to a stationary pattern observed numerically in m was proved in [23]. The 
objective of this work is to obtain asymptotic formulas for the pattern and the timing of its 
formation under the assumption of slow diffusion. 
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The model includes several time scales. The “spatial” variable is a scalar parameter x 
that measures the distance between the switching thresholds of a bacterium and characterizes 
its switching strategy (the smaller the x the more responsive is the strategy to variations of 
the environmental conditions). That is, a solution describes the evolution of the distribution 
of each of the two phenotypes over possible switching strategies parameterized by x. The 
diffusion of this distribution models sporadic changes of the switching strategy by individual 
organisms; it is assumed to be slow compared to the rate of growth and competition precesses. 
The rate of transition from one phenotype to the other in bacteria is assumed to be much 
higher than the rate of other processes. Every such transition is modeled by an instantaneous 
switch of a relay from state 1 to state —1 or vice versa (where each state represents a 
particular phenotype). 

The central assumption that we make is that the switching threshold for the transition 
of a relay (bacterium) from state 1 to state —1 is different from the switching threshold for 
the opposite transition (the difference of the thresholds 2x is called the bi-stability range). 
That is, we assume hysteresis, or multi-stability, in the switching response of bacteria to 
the exogenous stimulus. This hysteresis can be viewed as form of the persistent memory of 
past environmental conditions by organisms. Since the seminal work of Max Delbriik [ 13] . 
epigenetic differences arising in the process of cell differentiation have been attributed to 
multi-stability of living forms. A classical example of such multi-stability is the behavior 
of lac-operon in E. coli ( lac-operon is a collection of genes associated with transport and 
metabolism of lactose in the bacterium; expression of these genes can be turned on by 
molecules called inducers). Experimental studies of the regulation of enzymes in E. coli 
and yeast that were conducted as early as in the 1950-60’s effectively demonstrated that 
two phenotypes each associated with “on” and “off” state of lac-operon expression can be 
obtained from the same culture. The fraction of each phenotype in the total population 
depended on the history of exposure to the inducer, and each phenotype remained stable 
through multiple generations of the bacteria [5lll0[lll[l37ll30l[3Tll5l1J . This behavior resembles 
the two-threshold hysteretic relay shown in Fig. 12.11 Later, multi-stable gene expression 
and hysteresis have been well-documented in many natural as well as artificially constructed 
systems [4],[T5(!l9l|3T[[35] . In particular, recent experiments using molecular biology methods 
permitted a further study of the region of bi-stability of the lac-operon when multiple input 
variables are used to switch the lac-operon genes on and off. 

A substantial amount of experimental and theoretical work provides an evidence that 
organisms use various (often randomized) strategies of phenotype switching in order to adapt 
to changing environmental conditions (see [36] and references therein). It is intuitively clear 
that phenotypic diversity within the population can help to increase the chances of survival 
in varying environment. Switching strategies that establish phenotypic diversity are called 
bet-hedging [3], Indeed, experimental and theoretical models of adaptation demonstrate 
that bet-hedging can evolve to maximize the net growth of the total population HU2U2S1 
[3811321138] . In [IS], a simple differential model was used to show that switching strategies 
with two thresholds where the switching moment depends not only on the state of the 
environment, but also on the phenotype itself, can further increase Darwinian fitness of 


2 



species. These strategies exploit hysteresis and are described by a two-threshold relay shown 
in Fig. 12.II in the adiabatic limit when switching is faster than the characteristic rate of the 
environment variations. It was also shown in [18] that the optimal bi-stability range, that is 
the optimal separation of thresholds 2x of the relay, which maximizes the net growth rate 
of the population is different for different environmental inputs. For this reason, the model 
that is considered in this work includes the distribution of bacteria (relays) over a set of 
switching strategies which have different separation of switching thresholds with x taking 
the values from an interval [x, x \. The diffusion process included in the model acts as a factor 
that diversifies the switching thresholds (strategies) in the population, while the competition 
provides a mechanism of selection that may favor a subpopulation with a specific x for a 
given law of variation of the state of environment. More detailed discussion of the model 
assumptions and the biological background can be found in [17 ; . il8li22] . 

Parabolic equations with distributed non-ideal relays were previously used for modeling 
spatial structures that were observed in spatially distributed colonies of bacteria, for example 
concentric rings in Petri dish experiments p i24Il2TI[39] . In these models, all the relays have 
the same fixed separation of thresholds. The mechanism of pattern formation discussed 
in this work is based on the long-term memory of the hysteretic relays and is related to 
dynamics of free boundary separating the domain where relays are in state 1 from the 
domain where relays are in state —1 on the interval [x, x]. This mechanism is different from 
the mechanism found in [9[|26][27]. Other important aspects of dynamics of parabolic and 
hyperbolic differential equations with distributed relays were studied in p [71fl2lll6l[20ll2 1 11291 
l3mi331!491l5U] in relation to multiple applied problems (not related to patterns) such as phase 
transitions, population dynamics, multi-phase flows in porous media, thermostat control and 
others. 

In the following Sections 2 and 3 we present the model and a statement on its well- 
posedness that has been obtained earlier. Section 4 contains the main statement on the 
asymptotics of the patterns and its proof. 

2 Model description 

In this paper, we consider a class of models, which attempt to account for a number of 
phenomena, namely (a) switching of bacteria between two phenotypes (states) in response 
to variations of environmental conditions; (b) hysteretic switching strategy (switching rules) 
associated with bi-stability of phenotype states; (c) heterogeneity of the population in the 
form of a distribution of switching thresholds; (d) bet-hedging in the form of diffusion be¬ 
tween subpopulations characterized by different bi-stability ranges; and, (e) competition for 
nutrients. The resulting model is a reaction-diffusion system including, as reaction terms, 
discontinuous hysteresis relay operators and the integral of those. This integral can be in¬ 
terpreted as the Preisach operator [SI HB(131) 132],:34 ]|35]I33]III1 HSUSQj with a time dependent 
density (the density is a component of the solution describing the varying distribution of 
bacteria). In [23], we have shown that fitness, competition and diffusion can act together 
to select a nontrivial distribution of phenotypes (states) over the population of thresholds. 
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The main goal of this paper is to give a quantative description of this distribution for small 
diffusion. 

We assume that each of the two phenotypes, denoted by 1 and —1, consumes a different 
type of nutrient (for example, one consumes lactose and the other glucose). The amount of 
nutrient available for phenotype i at the moment t is denoted by fi(t) where i = ±1. The 
model is based on the following assumptions (see [22] for further discussion). 

1. Each bacterium changes phenotype in response to the variations of the variable w = 
/i/(/i + f- i) — 1/2, which measures the deviation of the relative concentration of the 
first nutrient from the value 1/2 in the mixture of the two nutrients. 

2. The input w = w(t) is mapped to the (binary) phenotype (state) of a bacterium 


r(t) = lZ x (w)(t), where 1Z X is the non-ideal relay operator with symmetric switching 


thresholds x, y = — x 
Section [3] 

with x > 0; 

i 

1 

see Fig. 12.11 and the rigorous definition (| 3 . 1 |) in 

.output 

' 
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x input 
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Figure 2.1: Non-ideal relay. 


3. The population includes bacteria with different bi-stability ranges (—x,x), where the 
threshold value x varies over an interval [x, x] C (0,1/2). We will denote by u(x, t ) the 
density of the biomass of bacteria with given switching thresholds ±x at a moment t. 

4. There is a diffusion process acting on the density u. 

5. At any particular time moment t, for any given x, all the bacteria with the switching 
threshold values ±x are in the same state (phenotype). That is, u(x,t ) is the total 
density of bacteria with the threshold x at the moment t and they are all in the same 
state. This means that when a bacterium with a threshold x' sporadically changes its 
threshold to a different value x, it simultaneously copies the state from other bacteria 
which have the threshold x. In particular, this may require a bacterium to change the 
state when its threshold changes. 

With these assumptions, we obtain the following model of the evolution of bacteria and 
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nutrients, 


Ut = Du xx + ^(1 + n x (w)) ufi + ^(1 - n x (w )) uf_ u 
< (1 + K x (w))udx, ( 2 . 1 ) 

f-i = ~f-i [ (1 -K x (w))udx, 

< X 

where u t and u xx are the derivatives of the population density u; D > 0 is the diffusion 
coefficient; dot denotes the derivative with respect to time; and all the non-ideal relays 7 Z x , 
x G [x, x], have the same input w — /i/(/i + /_i) — 1/2. We additionally assume the growth 
rate |(1 + ilZ x (w)) ufi based on the mass action law for bacteria in the phenotype i = ±1. 
The rate of the consumption of nutrient in the equation for f t = fi(t) is proportional to the 
total biomass of bacteria in the phenotype i, hence to the integral; x and x are the lower 
and upper bounds on available threshold values, respectively. 

We assume that a certain amount of nutrients is available at the initial moment; the 
nutrients are not supplied after that moment. We assume the Neumann boundary conditions 
for u, that is no flux of the population density u through the lower and upper bounds of 
available threshold values. 


3 Rigorous setting of a well-posed model 

3.1 Rigorous setting 

Throughout the paper, we assume that x G [x, x\ C (0,1/2). 

We begin with a rigorous definition of the hysteresis operator 1Z X ( non-ideal relay ) with 
fixed thresholds ±x. This operator takes a continuous function w = w(t) defined on an 
interval t G [0, T) to the binary function r = lZ x {w) of time defined on the same interval, 
which is given by 


K*W(t) = { 


( -1 
1 

r 0 


if w(t) < —x for some r G [0, t] 
and w(s) < x for all s G [r, t], 
if w(t) > x for some r G [0, t] 
and w(s) > — x for all s G [r, t\, 
if —x < w(t) < x for all r G [0, t\, 


(3.1) 


where is either 1 or —1 (initial state of the non-ideal relay 1Z X ). Since ro may take different 
values for different x, we write r 0 = r 0 (x). The function r 0 = r 0 (x) of x G [x, x\ taking values 
±1 is called the initial configuration of the non-ideal relays. In what follows, we do not 
explicitly indicate the dependence of the operator 1Z X on r 0 (x). 

In this paper, we assume that ro(x) is simple, which means the following. There is a 
partition x = x 0 < X\ < ■ ■ ■ < xjv 0 = x of the interval [x, x] such that the function r 0 (x), 
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which satisfies |r 0 (a;)| = 1 for all x < x < x, is constant on each interval (x k -i,x k ] and has 
different signs on any two adjacent intervals: 


(3.2) 


ro(x) = r 0 (x k ), x e [x k _ u x k ], k — 1,..., N 0 , 
r 0 (x fc -i)r 0 (x fc ) = -1, k = 2,..., N 0 , 

where the second relation holds if Nq > 2. 

We define the distributed relay operator TZ(w) taking functions w = w(t) to functions 
r = r(x, t ) by 

r(x,t) — lZ(w)(x,t) \=lZ x [w){t). (3.3) 

The function r(-,t) will be referred to as the configuration (state) of the distributed relay 
operator at the moment t. 

We set 


U(t) = / u(x,t)dx, V(u,w)(t) — / u(x, t)TZ x (w)(t) dx. 


(3.4) 


Here the first integral is the total mass of bacteria; V is the so-called Preisach operator with 
the time dependent density function u. Further, we replace the unknown functions j\ and 
/_i in system (12.11) by the functions v — /) + f_i (total mass of the two nutrients) and 
w = /i/(/i + /_i) — 1/2 (deviation of the relative concentration of the first nutrient from 
the value 1/2). Substituting 


fi = ( 2 + w ) v ' 


f -1 = ( 2 “ w ) v 


into the first equation of (12.ip . we obtain the equation 

T 


u t = Du x 


+ wTl(w) uv 


Furthermore, summing the second and the third equations of system (12.1|) and using the 
relationships v — j\ + /_ i, 2 wv — j\ — /_ i and the notation (13.41) . we get 


v = 


U „ / 4 1 

— + wV(u, w) ) v. 


Finally, from the second and the third equations of (12.11) . it follows that 


w = 


fif-i ~ fif-i _ fif-i 

(A + /-1) 2 “ 


uTZ x iw) dx = — 


+”) Q-") 


V(u, w) 


Thus, the resulting system, which is equivalent to equations (12.ip . takes the form 

T 


u t = Du x 


+ wlZ(w) uv, 


v = 


w = 


+ wV(u, w ) ) v, 


(3.5) 


+ ”) \l- w ) 


V(u,w), 
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where we assume the Neumann boundary conditions 


Un 


r—r. x=x 0 


(3.6) 


and the initial conditions 

u(x, 0) = u(0) = Uq, w(0) 


pX 

u 0 (x) >0, / uo(x) dx = 1, v 0 > 0, 

J X 


= w 0 , r(x, 0) = r 0 (x), 

|Wq | < x, r(x) is simple. 


(3.7) 


3.2 Well-posedness 

The problem (I3.5I) - (I3.7I) . which contains a discontinuous distributed relay operator, was 
shown in [22j to be well posed. We briefly summarize this result before proceeding with the 
analysis of long time behavior. 

Set Qt = (x,x) x (0,T) for T > 0. We will use the standard Lebesgue spaces L 2 (Q T ) 
and L 2 = L 2 (x,x)] the Sobolev spaces W 2 = W%(x,x), fceN; the anisotropic Sobolev space 
W^iQr) with the norm 

/ t T \ 1/2 

u(-,t)\\ 2 W 2 dt + J \\u t (-,t)\\l 2 dt\ ; 

o / 

and the space 

W(Qt) = W*'\Q t ) x C^T] x C'[ 0,T]. 

Definition 3.1. Assume that (woWo^o) £ W) x We say that ( u,v,w ) is a ( global) 
solution to problem (13.5h (13.7(1 if, for any T > 0, ( u,v,w ) £ W(Qr), r(-,t) is a continuous 
./^-valued function for t > 0, and relations (13.511 (13.7(1 hold in the corresponding function 
spaces. 

The following result was proved in [22] , 

Theorem 3.1. If (u 0 ,v 0 ,w 0 ) £ W\ x M 2 , then 

1. Problem (13.511 - (13.711 has a unique solution ( u,v,w ); 

2. The state r(-, t ) = IZ(w)(-, t ) of the distributed relay operator is simple for all t > 0; 

3. We have 

U_(t) > 0, (7(f) —* 1 + v 0 , 0 < u(-jt) —> ° in C[x, x], (3.8) 

x — x 

v(t) < 0, v(t)<v oe _Mt with p = 1/2 —x>0, (3.9) 

\w(t)\ <x< 1/2, (3.10) 

where the convergence takes place as t —> oo. 


m 


W a a,1 «T) 
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The behavior given by (13.8ft (13.101) is to be expected. Indeed, as we assume no supply of 
nutrients after the initial moment, the total amount of nutrients v(t) converges to zero. When 
the density of nutrients vanishes as a result of consumption by bacteria, the equation for the 
density u approaches the homogeneous heat equation with zero flux boundary conditions, 
which explains why the density of bacteria u(x,t) converges to a uniform distribution over 
the interval [x, x] as a result of the diffusion. 

Below we are interested in the limiting behavior of the configuration function r(x, t) given 
by (13.3p . or, in other words, in the limiting distribution of phenotypes over the population 
of thresholds. In [23] . we have proved that r(x,t ) converges to a step like profile as t —> oo, 
see Fig. 15.11 In the next section, we give a quantative description of this phenomenon under 
the assumptions that the initial amount of nutrient Vo and the diffusion coefficient D tend 
to zero, while the initial density Uq{;x) tends to the delta function. 

4 Fronts asymptotics 

4.1 Terminology 

A point Xj = Xj{t) which separates an interval on which the relays are in state 1 from an 
interval on which the relays are in state —1 is called a front ; cf. (13.21) . The total number of 
fronts can vary, but, by Theorem 13.11 remains finite at all times (equivalently, the state of 
the distributed relay operator remains simple at all times). A front can either stay (a steady 
front) or move right. That is, any Xj(t) is a non-decreasing function on the time interval of 
its existence. A front disappears if it hits another front or is hit by another front. Assume a 
front Xj(t) exists at a moment to- It is called immortal up to time T, T > to, if it does not 
disappear on the time interval to <t <T. 

Our goal is to give asymptotic formulas for the time moments at which steady immortal 
(up to time T) fronts appear as well as for their positions. This will be done as D —> 0, 
where D is the diffusion coefficient in (13.5ft responsible for random fluctuations of hysteresis 
thresholds of individual bacteria. 

4.2 Assumptions 

From now on, we fix a time interval [0,T], T > 0. Along with conditions (13.7ft . we 
that the initial density satisfies the relation 

Uo(x) < e, x G (x, x — e), 

where e > 0 will be chosen small enough and depending on D. 

We introduce the function 

F(x) 1/4 — (x — x) 2 

and set 

*+ £ dx 
F{xY 



assume 

(4.1) 

(4.2) 

(4.3) 
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(4.4) 


Set 


F := 



dx 

F( x y 


rx 

U(x,t) : = / u(y,t) dy , 

J X 


E(y):= 




dz. 


Using [2S Lemmas 4.1 and 5.4] and the fact that 0 < U_(t) — 1 < vq (see (j3.7[) and (I3.8|) ). 
we can choose positive functions e(D) and v 0 (D) such that e(D),v 0 (D ) —> 0 as D —» 0 and 


sup 

i£[ s l/2/2,T] 


U(-,t) - E ( 


^(Di) 1 / 2 ) 


= o(l) as D —* 0, 


C\x,x] 


(4.5) 


provided that uq{x) satisfies (14. ip with e = e(D) and that vq = vo(D). 

The main assumptions will be as follows. 

1. The initial density u 0 (x ) satisfies (14.ip with e = e(D), where e(D) is the above function. 

2. The initial amount of food is small: v 0 = vq(D), where Vq(D) is the above function. 

3. The initial value of the input is close to x: wq = wq{D) 1 where wq(D) is an arbitrary 
function such that wq{D) < x and wq{D) —» x as D —> 0. 

4. The initial conhguration is r 0 (x) = 1. 

Our goal is to determine the consecutive time moments t n at which the moving fronts 
become steady and their positions x n (n — 1,2,...) at these moments. 

We begin with the following recursive algorithm for determining (finitely or infinitely 
many) positive numbers 


0 Sq Si S 2 s n . . . , 2 /l, 2/25 ■ • • j Urn ■ • • 


(4.6) 


4.3 Algorithm 

Basis. Set Si := F, and let y\ be the (unique) root of the equation 2 E(y) = 1, y G (0, oo). 

Set G 2 (t) := —2 E (^y\ ^ • Note that G 2 (si) + 1 = 0 and G 2 (t) > 0 for all t > 0. 

Inductive conjecture. Fix n > 2. Assume that we have defined the sequences 
Si, ..., s n _i and yi,, y n -\ such that the function 

1 / i/qX 

G„(():=2g(-l)"«S^(|) ) (4.7) 

satisfies 

G n (s n -i) + 1 = 0, G n (t ) + 1 > 0, t > s n _i. (4-8) 

As we have seen, relations (14. 8 p hold for n — 2. 
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Below we will also use the equivalent recursive representation of G n (t ): 


Gift) = 0, G n (i) = -G„_,(()-2Eh„_ 1 (^) 1/2 V n> 2. 


(4,9) 


Inductive step. Now we will determine s n , y n , and the function G n+ i(t) satisfying 
relations (14. 9 p with n replaced by n + 1. 

1. Consider the equation 


(G n (t) + l)dt — F, s G [s n _i, oo). 


(4,10) 


' &n — 1 


Due to (14. 8p and the fact that G n (t) + 1 —* 1 as t —* oo, the left-hand side is a strictly 
increasing function of s that tends to oo as s —* oo. Therefore, Eq. (14.10p has a unique 
root, which we denote by s n . 


2. Consider the function 


H n (t, z ) := - G n (t ) - 2 E ^-jj/ 2 ) + 1 > f z - 0 - 


Relations (14.81) imply that 


On the other hand, due to 


— G n (t ) + 1 < 2, t > s r 


(4.11) 


(4.12) 


—G n (t) + 1 — G n _i(t) + 2E j ^ + 1, t > s n , 

and, hence, 

— G n (t ) + 1 > 0, t > s n , (4-13) 

due to the inductive conjecture (see (14.8p with n replaced by n — 1). 

Thus, we see that the equality H n (t, z) = 0 uniquely defines a smooth positive function 
z — Z n (t), t > s n , such that 

H n (t,z)> 0, z < Z n (t)] H(t,Z n (t)) = 0; I7 n (t,^)<0, z > Z n (t). (4.14) 


3. Now we have several possibilities depending on the geometry of the curve Z n (t) defined 
by (J4H: 

(a) Z n (t) > 0 for t > s n . In this case, we set s n := s n . 

(b) There is s n > s n such that Z n (t) < 0 for t G (s n , s n ) and Z n (t) > 0 for t > s n . 

(c) Otherwise, we terminate the process at the finite sequences {si,..., s n _i} and 

{l/l; ■ ■ ■ j Un— l}- 
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In cases |3a] and I3hl we set 


(4.15) 


_ ^n(^n) 

y™ 1/2 ‘ 

Sn 

In these two cases, it remains to check that the function 

G n+1 (t) := —G n (t) - 2 E (y n (y) V ^ (4.16) 

satisfies relations (14 .8 [) with n replaced by n + 1. This will guarantee that we can repeat 
the inductive step of the algorithm from section EOl with n replaced by n + 1. 

Combining (j4. 16[) with (144 ip . (14.14j) . and (144 5p . we have 

C ri+ i(s n ) + 1 = — G n (s n ) — 2 E ( y n ) + 1 — H ( s n , Z (, s n )) = 0. 

On the other hand, in cases l3al and l3bl we have Z n (s n ) < Z n (t) for t > s n . Therefore, 
relations (144 ip . (14451) . and (I4.16P imply that 

G n +\(t) + 1 = H n (t , Z n (s n )) >0, t > s n . 


4.4 Main result 

4.4.1 Formulation of the main result 

Assume that, for some N > 1, we have the sequences {si,... ,sv} and {y i,... ,yv} con¬ 
structed according to the above algorithm, and assume that sn < T, where T was fixed in 
Sec. 1441 

Theorem 4.1. For all n = 1,..., N, the n-th moving front becomes steady and immortal 
(up to time s^) at a moment 

t n = s n + o( 1) (4.17) 

and its position at this moment is 

1 /9 

x n = x- q n , where q n = 2 (y n + o(l)) (Ds n ) ; (448) 

here o(l) stands for functions of D that tend to 0 as D —> 0. Furthermore, x n < x^-i < 

■ ■ ■ < X\ < x. 


4.4.2 Preliminary discussion 

In section [5j we will prove Theorem 1441 Assume we have constructed n —1 fronts that became 
steady at the moments ti,, t n -1 at the positions x ±,..., x n -\ given by (14.17ft and (I4.18p . 
respectively. We will consider the formation of the n-th front on the time interval [t n -i,t n . 
We consider the case of an even n > 2. The case of odd n is analogous; see also Remark [54] 
below concerning n = 1. 
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Below we will consecutively consider time moments 

[tn-l — Sn-1 + 0(1)] < fn -1 < t n - 1/2 

< [f n = S n + 0(1)] < [r n = S n + 0(1)] 

< [On = S n + 0(1)] < [t n = S n + 0(1)] < r n 

which are characterized by the following properties: 

1. The function w(t) is increasing and the relays do not switch during the time interval 
t E (t„—i,t n — 1 / 2 ); furthermore, 


w(t n - 1 ) = -x + q n - 1 , w(t n - 1 / 2 ) = x. (4.19) 

2. The function w(t) is increasing and some relays switch during the time interval (£„_i/ 2 ,f n ); 
further, for some Z*(D) = o(l) (which will be chosen below to satisfy relations (I5.15P ). 
we have 

w(f n ) = x — 2Z*(D). (4.20) 

Note that, although the distance between w(f n ) and x is already 2 Z*(D) = o(l), it is 
still much larger than the desired distance 2 (y n + o(l)) (Dsn) 1 /" due to (15.121) . Also 
the moment f n , which is of order s n , may be much smaller than the desired moment, 
which is of order s n . 

3. The function w(t) may oscillate and the relays may, though not necessarily, switch 
during the time interval (r n ,r n ). The moment r n = s n + o(l) is already of the correct 
order, but the value w{r n ) may still have a wrong asymptotics. However, it satisfies 

w{f n ) < w(r n ) < x -2 (y n + o(l)) ( Ds n ) l/2 , (4.21) 

i.e., it is bounded from above by the correct asymptotics. 

4. The function w(t) may oscillate and the relays may, though not necessarily, switch 

during the time interval ( r n ,O n ). The moment O n = s n + o(l) is of the same (correct) 

order as r n , but now, additionally, the value w(O n ) has the correct asymptotics 

w(O n ) = x-2 (y n + o(l)) (Ds n ) 1/2 . (4.22) 

5. The function w(t) may oscillate and the relays may, though not necessarily, switch 

during the time interval ( O n ,t n ). The moment t n = s n + o(l) is of the same (correct) 

order as r n and 9 n , and the function w(t) does not oscillate any more after the moment 
t n . More precisely, w(t) is decreasing and the relays do not switch during the time 
interval ( t n ,r n \. Moreover, the n-tli immortal front is formed at the moment t n in the 
sense that 

w(t n ) = x - q n = x n , w(t n ) = 0. (4.23) 
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5 Proof of the main result 


In this section, we will prove Theorem 14.II Recall that, without loss of generality, we assume 
that n is even. In sections I5.1H5.5I below, we will consider in detail the respective time 
intervals from items 1-5 of section 14.4.21 

5.1 Dynamics for w £ (— x-\- q n -\,x) as £ £ t n _i/ 2 ) 

We will see below that w(t ) is increasing for t > t n -\ and achieves the value w{t ) = x. Let 
Li- 1/2 be the moment when this happens: u>(t n -i/ 2 ) = x. Note that if w(t) G (—x + q n -i,x ) 
and w(t ) > 0, then the relays do not switch. 

Assume we have proved at the previous step (for n replaced by n — 1) that 

w(t ) >0 for t G (L l - 1 ,r„_ 1 ] (5.1) 

for some r n _i > s n _ 1 not depending on D (we shall do this for k = n in the end; see 115.331) ). 
As long as w(t) remains positive and w(t ) < x holds, the relays do not switch, see Fig l5.1l a. 
Hence, using (14. 5 p and (j4. 18[) . we can write the third equation in (13. 5 p as follows (we omit 
“u(x,t) dx ” for brevity): 



= (U(x,t) - 2U(x n _ 1 ,t) +2U(x n _ 2 ,t) - 2U(x 1 ,t)) ( 5 . 2 ) 

= (j - ^1 + 2 J^(-l YE [{ yj + o(l)) 1 ^ + o(l) j , t > t n _ 1 

(here and below, we consider t < T only.) 


X n -\ 

%2 

Xi 

42 

l 

IH 


X 

W(t) 

Xn -1 

x 2 

Xi 

X 



X 





X 


Figure 5.1: Spatial configuration of the relays lZ x (w) at a moment t > t n -\. a) The input 
satisfies w{t ) G (—x + q n -i, x). Hence, the relays do not switch and all the fronts are steady, 
b) The input satisfies w(t) G (x,x — q n ). Hence, the relays may switch and the n-th front 
W{t) can move. 

We set 

p(t) :=x + w(t). 
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Taking into account (15.2p . definition (14. 2 p of F(x), definition (14. 7 p of G n (t), and the fact 
that n is even, we see that p(t) satisfies 


p(t) = F(p)(G n (t) + 1 + fji(t,D)), t>t n - 1, (5.3) 

p(t n - 1 ) = 2(j/ n _i + o(l))(ds„_i) 1/2 , (5.4) 

where 

sup \n(t,D)\ < e n = e n (D) for a nonnegative e n (D) = o(l). (5.5) 

*£[ s l/2/ 2 T] 

Due to (14.81) and (15.51) . the right-hand side in (15.31) is positive for t G [r n _i,T], provided 
D is small enough. Together with (15. 11) . this means that the solution to (15.31) . (15.4[) increases 
until a moment f n _i/ 2 (if it exists) at which it achieves the value x + x. This proves that (j5.3[) 
is equivalent to the third equation in (13. 5 p not only for t G (t n _i,r n _i) but actually for 
t £ [fn— 1) Li— 1 / 2 )■ 

The following lemma determines the time moment £ n _ i/ 2 . 

Lemma 5.1. The equation 

rs px+x j 

/ {G n (t) + l)dt= ^ (5.6) 

Jsn-l JO r \ X > 

has a unique root on the interval s G (s n _i,oo). Denote it by s n _i/ 2 - Then t n _i/ 2 = 
S n - 1/2 + O(l) p(f n _l /2 ) > 0. 

Proof. Due to (14.81) and the fact that G n {t) + 1 —> 1 as t —» cx), Eq. (15.61) has a unique root 
s n - 1/2 • 

To prove the asymptotics for we denote by p±(^) the solutions of 

P±(t) = F(p±)(G n (t) + 1 ±e„), t > t n _ 1 , 

with the same initial data at t = f n _i as in (j5.4[) . Obviously, 

P-(t) <P( f ) <P+(t). (5.7) 

Let tn-i/ 2 ,± be the first moment at which p±(t) = x-\-x. Then, integrating the differential 
equation for p±(t), we obtain 

i/2,± rx+x j 

(G n (t) + l±e n )dt= / 

+ 0 ( 1 ) J 0 ( 1 ) x 

Using the implicit function theorem in a neighborhood of the point (i„_ 1 / 2 ,±j D) = (s„_i/ 2 , 0), 
we obtain t n -i/ 2 ,± = s n - 1/2 + o(l). Therefore, taking into account (I5.7[) . we have f n -i /2 = 
Sn—1/2 + o(l). 

The inequality p(£ n -i/ 2 ) > 0 has been proved before the lemma. □ 
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5.2 Dynamics for w £ (x, x — q n ) as t E (t n _ i/ 2 , f n ) 

By Lemma EU tu(in-i/ 2 ) > 0. Therefore, as long as ii;(t) remains positive, the input w(t) 
increases and hence switches the relays. However, if w(t) becomes negative, the relays will 
stop switching. To describe the dynamics of w(t), we introduce the function 


W(t) := max w(s), 
se [t n - 1/2 ,t] 


where t n _i /2 is defined in Lemma E~T1 see Fig IS.ll b. Then, similarly to (15.21) . we obtain from 
the third equation in (13. 5 ft and from relation (14.5ft 


w(t) = -[l~ w2 


" W(t ) rx n -1 rx „-2 rXn -3 

- + - +••• + 

c JW(t) JXn-1 J Xn-2 


'XI 


(i ~~ W ^j _ 0 + ( x n- l, t) ~ 2(7(x n _ 2 , t) H-(-2(7(xi, f)) 


- - w 2 1 I 1 — 2E 


Setting 


we have 


(|^)^E(-i)^(te + o(i))(f) 1/2 


q(t) := x — w(t), Q(t) min q(s) = x — W(t), 


+ o(l) , t > ; 


g(f) = F(g) -G n (f) - 2£ 


Q 


2(Df) 1 / 2 


+ 1 + D) ) , t > t n —\/2i 


where 


with S 1/2 defined in 


sup |/i(f, _D)| —* 0 as .D —> 0 

[si/ 2 /2,X'] 

We consider Eq. (15.8[) with the initial condition 
q(t n - 1 / 2 ) =x-x. 


(5.8) 

(5.9) 

(5.10) 


We begin with the following observation, which shows that the (n — l)-th front is immortal. 
Lemma 5.2. The solution q(t ) to problem (15.8ft . (15.101) cannot achieve the value q n ~i- 

Proof. Assume, by contradiction, that q(t) achieves the value g n _ 1 (= 2(?/ ri _ 1 +o(l))(.Ds n _i) 1 / 2 )) 
for the first time at some moment t = t* (> y 2 > s n _ 2 ). At this moment Q(t*) = q(t*) = 
q n - 1 - Since t* > s n _ 2 , it follows from (14.8ft (with n replaced by n — 1) that G n -i(t*) + 1 > 0. 
Therefore, using (15.81) and (14.9ft . we have 

q(t*) — F(0)(G n -i(t*) + l) + o(l) > 0 


for all sufficiently small D, which is impossible. 


□ 


n—l/2 • 
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To proceed, we will need a specific (positive) function Z*(D) such that 


Z*(D) —^ 0 as D —> 0, 

Z*(D) 

D 1 / 2 


—> oo as D —> 0. 


To define it, we introduce the new function z(t) : = 


Q(t) 

2D 1 / 2 


and set 


Z(t) := min z(s). 
se[in- 1/2A 


(5.11) 

(5.12) 

(5.13) 


Then we consider Eq. (15.8ft for t > s n , where s n is the time moment defined in the inductive 
step CD of the algorithm from section l4~3l For such t, it is equivalent to 


z(t) = 


F(0)(H n (t,Z) + o(l)) 
2D 1 ! 2 


(5.14) 


where H(t, •) is given by (14.111) . We choose the desired Z*(D) to satisfy, for all t G [s n ,T], 

(5.15) 


H n (t, C) + o(l) >0 for C < Z n (t) - Z*(D), 

H n {t, () + o(l) < 0 for C > Z n (t) + Z*(D), 

F(0)\H n (t, C) + o(l)| < 2h n Z*(D) for ( e [Z n (t) — Z*(D), Z n (t) + Z*(D)], 


where o(l) is the same function as in (15.141) and h n > 0 does not depend on t and D. This 
can be done because dH n (t, Z n {t))/dQ ^ 0 for all t G [S n , T], 

The following lemma defines a time moment f n satisfying (14.201) . 

Lemma 5.3. There is a time moment f n = s n + o(l) such that q(t ) < 0 for all t e [t n _ 1 / 2 , f n ] 
and q{j n ) = 2 Z*(D). 

Proof. As long as q(t) < 0 and Z„(D) < q(t) < x — x, we have Q(t) = q{t). Therefore, we 
can rewrite Eq. (15.8p as follows: 


q{t) = —F(q) ( G n (t ) + 1 + nx (t, q, D)) , t > t n — 1 / 2 , 
where, due to (15.9ft and (15.12ft . 


(5.16) 


sup sup |/ii(t, q, D)\ <S — 5(D) for a nonnegative 5(D) = o(l). (5.17) 

<je[0,x-x] tg[si/ 2 / 2 ,T] 

Due to (14.8P and (|5.17|). the right-hand side in (f5TT6|) is negative for t G \fn— 1/2 j T \, 
provided D is small enough. Hence, q(t) < 0 for t > t n _i/ 2 - 
Denote by q±(t) the solution to the equation 


q±(t) = —F(q±) ( G n (t) + 1 ±5), t> t n _ 1 / 2 , 
with the same initial condition as in (15. lOjl . q(t n - 1 / 2 ) = x — x. 


(5.18) 
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Then, we have <?_(£) < q{t ) < q+(t). Therefore, it suffices to show that q±(t) achieves the 
value 2D 1 / 2 ~ x for the first time at a moment f± = s n + o(l). Equation (j5.18j) yields 


/ (G n (t) + l±6)dt= —— 

J J2Z,{D) P K x 

s n- l/2+°(l) 


Using the fact that s n is a root of (14.101) and s n _i/ 2 is a root of (15.6[) as well as the 
definition of F in (14.41) and the symmetry of F(x), we have 


s n-l/2 


(G n (t) + 1 — / ( G n [t ) + 1 


s n- 1/2 




— F — 


dx 

Ffx) 


Sn — 1 
px—x 


(G n (t) + 
dx 


F(x) 


Hence, the implicit function theorem in a neighborhood of (r±, U) = (s n , 0) together with (15. lip 
yield f± = s n + o(l). □ 


5.3 Dynamics for w £ (x, x — q n ) as f 6 (f n , r n ) 

Assume that case [3b] in the inductive step of the algorithm from section 14.31 holds with some 
s n > s n . If case I5al holds, we omit this step and proceed with section I5~41 below. 

The following lemma defines a time moment r n satisfying (14.211) . 

Lemma 5.4. There is a time moment r n = s n + o(l) < s n such that q(r n ) < 0 and 

2(Z n (r n ) - Z^D))D^ 2 < q(r n ) = Q(r n ) < 2Z*(D). (5.19) 

Proof. To follow the proof, we refer the reader to Fig. 15.21 We consider equation (15.14j) with 
the initial data 

z{f n ) = Z if [D)D~ 1 G. (5.20) 

First, note that z(t) > Z(t ) > Z n (t) — Z*(D) for all t e [T n ,s n \. Indeed, if the equality 
Z(t) = Z n (t) — Z*(D) is achieved for the first time at some moment t G [f n ,s n ], then 
z(t) = Z(t) = Z n (t) — Z*(D). Moreover, since Z n (t) — Z*(D) is nonincreasing, we have 
z(t) < 0 at this moment. On the other hand, relations (15.14[) and (15.151) imply that z(t) > 0. 
Thus, the first inequality in (15.19[) holds with r n replaced by any t < s n . 

Furthermore, Z(t) is nonincreasing. Hence, Z(t) < Z(f n ) < Z*(D)D~ 1 t 2 for all t > f n 
by Lemma [5731 It remains to find r n such that z(r n ) = Z(r n ) and q(r n ) < 0. 

Next, we introduce a positive function a*(D) such that a*(D) —> 0 and 

Z n (t) < -2 y/ZjD), t<s n - a*(D). (5.21) 
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Figure 5.2: Dynamics for t close to s n . Three black curves are the graphs of Z n (t) and 
Z n (t) ± Z*. The upper grey curve is the graph of the solution z(t) to equation (15.14H . The 
lower (horizontal) grey curve is the graph of Z(t) given by (j5. 13f) . 


Consider the time moment 


r* = r*(D) := s n - a*(D) - \fzjj5). 


First, assume that 

z(t*) > Z(r*) (5.22) 

(the case z(r*) = Z(t *) will be considered at the end of the proof). Denote by r 0 the time 
moment preceding r* such that 

z(t 0 ) = Z (t 0 ), z(t) > Z (t) = Z (r 0 ), t e [to,t*]. (5.23) 

Obviously, Z(r 0 ) G [Z n (r 0 ) — Z*, Z n (r 0 ) + Z*]. Indeed, otherwise, Z(r 0 ) > Z n (r 0 ) + Z* and, 
due to (I5.15P and equation (15. 14jl . we would have z(r 0 ) < 0. 

Note that, due to (I5.2ip . we have Z(r 0 ) > Z n (r 0 + \fZ~*) + Z*. Denote by 


r oi £ [hu r o + v^] 


a moment such that 

Z(r 0 ) = Z n (r 0 1 ) + Z*. (5-24) 

If Z(r 0 1) < Z(r 0 ), then there exists a desired time moment r n G [t 0 ,t 01 ] and the proof is 
complete. 
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Assume that Z(r 01 ) = Z(r 0 ). Then (I5.24p yields 

Z(t 01 ) = Z n (T 01 ) + Z*. (5.25) 

Let us estimate the velocity z(t) for t E [t 0 ,t 01 ]. For such t, we have Z(t) > Z n (t ) — Z * 
and, due to (I5.24p . Z(t) = Z(r 0 ) < Z n (t ) + Z*. Hence, (I5.15P implies that z(t) < h n Z*D~ 1 / 2 
for t E [to,t 0 i]. Therefore, 

z(r 0 i) < Z(t q ) + h n Z*\fzlD~ x t 2 . (5.26) 

Set 

7"02 := Toi + y/Z: 

From Z(t ) < 0 it follows that Z(t) > Z n (t) + Z* for t E [r 0 i,r 0 2 ]. If Z(r 0 2 ) < Z(r 0 i), then 
there exists a desired time moment r n E [ 701 , 702 ] and the proof is complete. 

Assume that Z(r 0 2 ) = Z(roi). It follows from (15.21ft and (15.251) that 

Z{tq 2 ) = Z(roi) = Z n (7oi) + Z* > Z n (ro2) + Z* + 2Z*. (5.27) 

Furthermore, due to (15. 15p and equation (15.1 4ft . we have 

z{t) < 0, i £ [%,%]. 

Therefore (see (15.26ft ). 

z(r 02 ) < z(t 0 i) < Z(t 0 ) + h n Z* v / XL rl/2 . (5.28) 

Now inequality (15.27ft . relations (15.15ft . and equation (15.14ft imply that z(t) < —hnZ^D^ 1 ^ 2 
at least as long as z{t) > Z(r 02 ) = Z(r 0 ), where h n does not depend on t and D. Therefore, 
taking into account (15.28ft . we see that z(t) will achieve the value Z(r 0 2 ) = Z(r 0 ) at a time 
moment r n E [r 0 2 , 7- 02 + h n y/Zl/h n ]. If follows from (15.23P and from the dehnitions of r*, r 0 , 
tqi, r 0 2 , and r n that r n E (r*,s n ). Hence, r n is the desired time moment. 

Finally, if (I5.22p is not valid, then z(r*) = Z(r*) and, as we already know, Z(r*) < 
Z*/}” 1 / 2 . The equality z(t*) = Z(t *) implies q(rf) < 0. If g(r*) < 0, then we can choose 
T n = r*. If g(r*) = 0, then relations f!5.14p . (15.15P imply Z(r*) E [Z n (r*) — Z*, Z n (r*) + Z*] 
and all the argument following Eq. (15.23P can be repeated for To = r*. □ 

5.4 Dynamics for w £ (x, x — q n ) as t £ (r n ,0 n ) 

Lemma 5.5. There is a time moment 9 n = r n + o(l) such that q(t ) < 0 for t E [r n , d n ] and 
q(6 n ) = 2{y n + o(l)) (Dsn) 1 ^. 

Proof. 1. Setting z(i) = ; we consider the problem (see Lemma 15.41) 

f F(0)(ff(t,Z(f)) + o(l)) 

^ W 2D 1 /2 ’ (5.29) 

{ y n s]! 2 - Z,(£>) < z(r n ) < Z^ZT 1 / 2 , 
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+ ZfiD) + o(l), then 


where H(t,z ) is given by 04. lip . If z(r n ) < Z n {r n ) + ZfiD) = y n s l J 2 
we take := r n and complete the proof. 

Assume that z(r n ) > Z n (r n )-\-ZfiD). We need to prove that there exist positive functions 
e*(D) and 6*(D) that tend to zero as D —* 0 and such that the solution z(t) to 05.14p satisfies 

z(t) < 0, t G [r n ,T n + 4(D)], y n < z(r n + 4(D)) < (y n + e*(D))s* /2 . 

To construct the functions £*(D) and 4(D), we consider positive sequences £ k ,5 k —$■ 0 
such that 4 < s n / 2. Due to 04.1311 . 04.14p . and 04.15p . there exists a sequence c k = c k (£k ) 
such that 

D(t, z) < -Ck, t e [s n , T], z> Z n (t) + £, fc. (5.30) 


Then 05.29p and 05.301) imply that ^(t) achieves the value Z n {t) + £k at a time moment 6 n k 
that satisfies 

ZJD)D- 1/2 AZJD) 

0 <9 nk ~T n < \\ . . , < ^ 

“ inf \z(t) ~ F(0)c k ~ 

te[s,2s n ] 

for all D < Dk, where Dk is a strictly decreasing sequence with Z\. —> 0. Therefore, 

1 /9 

= z n (9 nk ) + = ( 2 / n + £fc)s n where -> 0 as k ->■ oo. Now we set £*(D) := i k and 

:= 4 for D G (D fc , □ 

5.5 Dynamics for w E (x,x — q n ) as t E ( 9 ni t n ) 

To describe the dynamics for t > 9 n , we consider 05.8p with the initial data 


q(9 n ) = 2{y n + o(l)) (Ds n ) 


1/2 


(5.31) 


(see Lemma l5.5p . The following lemma will allow us to justify 05. ip (with n — 1 replaced 
by n) and thus to complete the proof of Theorem 14.11 

Lemma 5.6. There exists 9 n = s n + o(l) > 9 n and r n (> s n ) not depending on D such that 
the solution q(t ) to 05. 8 1) . (15.3 ip satisfies q(t ) > 0 fort G [9 n ,r n ] and w(r n ) > —x. 

Proof. Since F(q) > 1/4 — x 2 for q G [0,x — x], it suffices to estimate the expression in 
brackets in 05.80 . Using the definition of Q(t), the monotonicity of E(-), and relations 05.3ip 
and 04.81) (with n — 1 replaced by n), we see that there exists 9 n — s n + o(l) > 9 n such that 


G n {t) - 2 E ( 2 (^ 172 ) + 1 + D) 

> - G n ( t ) - 2 E ^2^^172) + 1 + jD ) 

= G n+ i(t) + 1 + o(l) > 0 for t > 0 n . 


(5.32) 


Equation 05. 8 p and inequality 05.32p imply that q(t) > 0 for t > 9 n . Obviously, we can 
choose r n > s n independent of D such that w(r n ) =x — q(r n ) > —x. □ 
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Finally, we determine the time moment t n (see (14.2311 ) at which the n-th front becomes 
steady and immortal and its position x n = x — q n as follows: 

t n ■■= sup {t e [9 n , 9 n \ : q(t) = 0, q(t) = Q(t)}, q n := q(t n ) = Q(t n ). 

Due to Lemmas 15.51 arid 15.61 t n = s n + o(l) and 

q(t) >0 for t G (t n ,r n ], (5.33) 

This justihes (15.11) (with n — 1 replaced by n). Moreover, since q(t n ) = 0, the value y : = 

, satishes (due to (15.81) and (I5.9R ) 

2 (Dt n y/- 

0 = —G n {s n + o(l)) — 2 E(y) + 1 + o(l). 

On the other hand, the dehnition (14.111) of H n (t,z), the definition (j4.14[) of Z n (t) and the 
dehnition (j4.15[) of y n imply that y n is the root of the equation 

0 = — G n (s n ) — 2 E(y n ) + 1. 

Therefore, y = y n +o(l), which implies q(t n ) = Q(t n ) = 2(r/ n +o(l)) . This completes 

the proof of Theorem 14.1[ 

Remark 5.1. An essential ingredient in the proof is the approximation (j4.5j) of the integral 
of the solution U by the error function E(-). It is valid for t separated from zero and thus 
can be used for n > 2. In the case n = 1, the analog of Eq. (15.31) will be q(t) = f(q)U_(t) on 
the interval t G (0,fi/ 2 ), where q(t) =x — w(t), q( 0) = o(l), and q{ti/ 2 ) = x + x. On this 
time interval, one can use the approximation U_(t) = 1 + o(l) (see (13.71) and (13.81) ) instead 
of (g3D. 

6 Numerics 

Figure 16.11 illustrates the values of s n and q n for n = 1,..., 10 from Theorem 14.11 found 
numerically. We note that s n = s n for n = 1,..., 6 (i.e., case [3a] from the inductive step of 
the algorithm from section 1473] takes place) and s n < s n for n — 7, ..., 10 (i.e., case [3b] from 
the inductive step takes place). In general, it is an open question whether one of the two 
cases [3a] and [3b] takes place for all n, or both may fail for some n and hence the resulting 
sequences (I4.6[) are finite. 
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Figure 6.1: The two graphs indicate the values of a) s n and b) q n for n = 1,..., 10 from 
Theorem 14.11 The maximal and minimal values for admissible thresholds are x = 1/100 and 
x — 1/4, respectively. 
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